#include "ManagedLasoProvider.h"

/* Table of constant values */

static integer c__0 = 0;
static integer c__1 = 1;

 /* Subroutine */int SmartMathLibrary::LasoPack::ManagedLasoProvider::dnppla_(S_fp
   op, S_fp iovect, integer *n, integer *nperm, integer *nop, integer *nmval,
   doublereal *val, integer *nmvec, doublereal *vec, integer *nblock,
   doublereal *h__, doublereal *hv, doublereal *p, doublereal *q, doublereal
   *bound, doublereal *d__, doublereal *delta, logical *small, logical *raritz,
   doublereal *eps)
{
  /* System generated locals */
  integer val_dim1, val_offset, vec_dim1, vec_offset, h_dim1, h_offset, hv_dim1,
    hv_offset, p_dim1, p_offset, q_dim1, q_offset, i__1, i__2, i__3, i__4;
  doublereal d__1;

  /* Local variables */
  /* extern doublereal ddot_(integer *, doublereal *, integer *, doublereal *, 
  integer *);*/
  static doublereal hmin, hmax, temp;
  //extern doublereal dnrm2_(integer *, doublereal *, integer *);
  static integer i__, j, k, l, m;
  //extern /* Subroutine */ int dcopy_(integer *, doublereal *, integer *, 
  // doublereal *, integer *);
  static doublereal dzero[1];
  //extern /* Subroutine */ int daxpy_(integer *, doublereal *, doublereal *, 
  // integer *, doublereal *, integer *);
  static integer jj;
  //extern /* Subroutine */ int dlaeig_(integer *, integer *, integer *, 
  // integer *, doublereal *, doublereal *, integer *, doublereal *, 
  // doublereal *, doublereal *, doublereal *, doublereal *, 
  // doublereal *, doublereal *, doublereal *);
  static integer kk;
  //extern /* Subroutine */ int dlager_(integer *, integer *, integer *, 
  // doublereal *, doublereal *, doublereal *);

  /* THIS SUBROUTINE POST PROCESSES THE EIGENVECTORS.  BLOCK MATRIX */
  /* VECTOR PRODUCTS ARE USED TO MINIMIZED THE NUMBER OF CALLS TO OP. */
  /* IF RARITZ IS .TRUE.  A FINAL RAYLEIGH-RITZ PROCEDURE IS APPLIED */
  /* TO THE EIGENVECTORS. */

  /* Parameter adjustments */
  q_dim1 =  *n;
  q_offset = 1+q_dim1 * 1;
  q -= q_offset;
  p_dim1 =  *n;
  p_offset = 1+p_dim1 * 1;
  p -= p_offset;
  hv_dim1 =  *nperm;
  hv_offset = 1+hv_dim1 * 1;
  hv -= hv_offset;
  h_dim1 =  *nperm;
  h_offset = 1+h_dim1 * 1;
  h__ -= h_offset;
  val_dim1 =  *nmval;
  val_offset = 1+val_dim1 * 1;
  val -= val_offset;
  vec_dim1 =  *nmvec;
  vec_offset = 1+vec_dim1 * 1;
  vec -= vec_offset;
  --bound;
  --d__;

  /* Function Body */
  dzero[0] = 0.;
  if (!(*raritz))
  {
    goto L190;
  }

  /* ------------------------------------------------------------------ */

  /* THIS CONSTRUCTS H=Q*AQ, WHERE THE COLUMNS OF Q ARE THE */
  /* APPROXIMATE EIGENVECTORS.  TEMP = -1 IS USED WHEN SMALL IS */
  /* FALSE TO AVOID HAVING TO RESORT THE EIGENVALUES AND EIGENVECTORS */
  /* COMPUTED BY DLAEIG. */

  i__1 =  *nperm **nperm;
  dcopy_(&i__1, dzero, &c__0, &h__[h_offset], &c__1);
  temp =  - 1.;
  if (*small)
  {
    temp = 1.;
  }
  m =  *nperm %  *nblock;
  if (m == 0)
  {
    goto L40;
  }
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    dcopy_(n, &vec[i__ *vec_dim1 + 1], &c__1, &p[i__ *p_dim1 + 1], &c__1);
    /* L10: */
  }
  (*iovect)(n, &m, &p[p_offset], &m, &c__0);
  (*op)(n, &m, &p[p_offset], &q[q_offset]);
  ++(*nop);
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    i__2 =  *nperm;
    for (j = i__; j <= i__2; ++j)
    {
      jj = j - i__ + 1;
      h__[jj + i__ * h_dim1] = temp * ddot_(n, &vec[j *vec_dim1 + 1], &c__1,
        &q[i__ *q_dim1 + 1], &c__1);
      /* L20: */
    }
    /* L30: */
  }
  if (*nperm <  *nblock)
  {
    goto L90;
  }
  L40: m +=  *nblock;
  i__1 =  *nperm;
  i__2 =  *nblock;
  for (i__ = m; i__2 < 0 ? i__ >= i__1: i__ <= i__1; i__ += i__2)
  {
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ -  *nblock + j;
      dcopy_(n, &vec[l *vec_dim1 + 1], &c__1, &p[j *p_dim1 + 1], &c__1);
      /* L50: */
    }
    (*iovect)(n, nblock, &p[p_offset], &i__, &c__0);
    (*op)(n, nblock, &p[p_offset], &q[q_offset]);
    ++(*nop);
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ -  *nblock + j;
      i__4 =  *nperm;
      for (k = l; k <= i__4; ++k)
      {
        kk = k - l + 1;
        h__[kk + l * h_dim1] = temp * ddot_(n, &vec[k *vec_dim1 + 1],  &c__1,
          &q[j *q_dim1 + 1], &c__1);
        /* L60: */
      }
      /* L70: */
    }
    /* L80: */
  }

  /* THIS COMPUTES THE SPECTRAL DECOMPOSITION OF H. */

  L90: hmin = h__[h_dim1 + 1];
  hmax = h__[h_dim1 + 1];
  dlager_(nperm, nperm, &c__1, &h__[h_offset], &hmin, &hmax);
  dlaeig_(nperm, nperm, &c__1, nperm, &h__[h_offset], &val[val_offset], nperm,
    &hv[hv_offset], &bound[1], &p[p_offset], &d__[1], &q[q_offset], eps, &hmin,
    &hmax);

  /* THIS COMPUTES THE RITZ VECTORS--THE COLUMNS OF */
  /* Y = QS WHERE S IS THE MATRIX OF EIGENVECTORS OF H. */

  i__2 =  *nperm;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    dcopy_(n, dzero, &c__0, &vec[i__ *vec_dim1 + 1], &c__1);
    /* L120: */
  }
  m =  *nperm %  *nblock;
  if (m == 0)
  {
    goto L150;
  }
  (*iovect)(n, &m, &p[p_offset], &m, &c__1);
  i__2 = m;
  for (i__ = 1; i__ <= i__2; ++i__)
  {
    i__1 =  *nperm;
    for (j = 1; j <= i__1; ++j)
    {
      daxpy_(n, &hv[i__ + j * hv_dim1], &p[i__ *p_dim1 + 1], &c__1, &vec[j
        *vec_dim1 + 1], &c__1);
      /* L130: */
    }
    /* L140: */
  }
  if (*nperm <  *nblock)
  {
    goto L190;
  }
  L150: m +=  *nblock;
  i__2 =  *nperm;
  i__1 =  *nblock;
  for (i__ = m; i__1 < 0 ? i__ >= i__2: i__ <= i__2; i__ += i__1)
  {
    (*iovect)(n, nblock, &p[p_offset], &i__, &c__1);
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ -  *nblock + j;
      i__4 =  *nperm;
      for (k = 1; k <= i__4; ++k)
      {
        daxpy_(n, &hv[l + k * hv_dim1], &p[j *p_dim1 + 1], &c__1, &vec[k
          *vec_dim1 + 1], &c__1);
        /* L160: */
      }
      /* L170: */
    }
    /* L180: */
  }

  /* ------------------------------------------------------------------ */

  /* THIS SECTION COMPUTES THE RAYLEIGH QUOTIENTS (IN VAL(*,1)) */
  /* AND RESIDUAL NORMS (IN VAL(*,2)) OF THE EIGENVECTORS. */

  L190: if (!(*small))
  {
    *delta =  - (*delta);
  }
  m =  *nperm %  *nblock;
  if (m == 0)
  {
    goto L220;
  }
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    dcopy_(n, &vec[i__ *vec_dim1 + 1], &c__1, &p[i__ *p_dim1 + 1], &c__1);
    /* L200: */
  }
  (*op)(n, &m, &p[p_offset], &q[q_offset]);
  ++(*nop);
  i__1 = m;
  for (i__ = 1; i__ <= i__1; ++i__)
  {
    val[i__ + val_dim1] = ddot_(n, &p[i__ *p_dim1 + 1], &c__1, &q[i__ *q_dim1 +
      1], &c__1);
    d__1 =  - val[i__ + val_dim1];
    daxpy_(n, &d__1, &p[i__ *p_dim1 + 1], &c__1, &q[i__ *q_dim1 + 1], &c__1);
    val[i__ + (val_dim1 << 1)] = dnrm2_(n, &q[i__ *q_dim1 + 1], &c__1);
    /* L210: */
  }
  if (*nperm <  *nblock)
  {
    goto L260;
  }
  L220: ++m;
  i__1 =  *nperm;
  i__2 =  *nblock;
  for (i__ = m; i__2 < 0 ? i__ >= i__1: i__ <= i__1; i__ += i__2)
  {
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ - 1+j;
      dcopy_(n, &vec[l *vec_dim1 + 1], &c__1, &p[j *p_dim1 + 1], &c__1);
      /* L230: */
    }
    (*op)(n, nblock, &p[p_offset], &q[q_offset]);
    ++(*nop);
    i__3 =  *nblock;
    for (j = 1; j <= i__3; ++j)
    {
      l = i__ - 1+j;
      val[l + val_dim1] = ddot_(n, &p[j *p_dim1 + 1], &c__1, &q[j *q_dim1 + 1],
        &c__1);
      d__1 =  - val[l + val_dim1];
      daxpy_(n, &d__1, &p[j *p_dim1 + 1], &c__1, &q[j *q_dim1 + 1], &c__1);
      val[l + (val_dim1 << 1)] = dnrm2_(n, &q[j *q_dim1 + 1], &c__1);
      /* L240: */
    }
    /* L250: */
  }

  /* THIS COMPUTES THE ACCURACY ESTIMATES.  FOR CONSISTENCY WITH DILASO */
  /* A DO LOOP IS NOT USED. */

  L260: i__ = 0;
  L270: ++i__;
  if (i__ >  *nperm)
  {
    return 0;
  }
  temp =  *delta - val[i__ + val_dim1];
  if (!(*small))
  {
    temp =  - temp;
  }
  val[i__ + (val_dim1 << 2)] = 0.;
  if (temp > 0.)
  {
    val[i__ + (val_dim1 << 2)] = val[i__ + (val_dim1 << 1)] / temp;
  }
  val[i__ + val_dim1 * 3] = val[i__ + (val_dim1 << 2)] *val[i__ + (val_dim1 <<
    1)];
  goto L270;

} /* dnppla_ */
